home *** CD-ROM | disk | FTP | other *** search
- /* Floating-point printing for `printf'.
- This is an implementation of a restricted form of the `Dragon4'
- algorithm described in "How to Print Floating-Point Numbers Accurately",
- by Guy L. Steele, Jr. and Jon L. White, presented at the ACM SIGPLAN '90
- Conference on Programming Language Design and Implementation.
-
- Copyright (C) 1992, 1993, 1994 Free Software Foundation, Inc.
- This file is part of the GNU C Library.
-
- The GNU C Library is free software; you can redistribute it and/or
- modify it under the terms of the GNU Library General Public License as
- published by the Free Software Foundation; either version 2 of the
- License, or (at your option) any later version.
-
- The GNU C Library is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- Library General Public License for more details.
-
- You should have received a copy of the GNU Library General Public
- License along with the GNU C Library; see the file COPYING.LIB. If
- not, write to the Free Software Foundation, Inc., 675 Mass Ave,
- Cambridge, MA 02139, USA. */
-
- #include <ansidecl.h>
- #include <ctype.h>
- #include <stdio.h>
- #include <float.h>
- #include <limits.h>
- #include <math.h>
- #include <stdarg.h>
- #include <stdlib.h>
- #include <localeinfo.h>
-
- #include <printf.h>
-
- #define NDEBUG
- #include <assert.h>
-
- #define outchar(x) \
- do \
- { \
- register CONST int outc = (x); \
- if (putc (outc, s) == EOF) \
- return -1; \
- else \
- ++done; \
- } while (0)
-
- #if FLT_RADIX != 2
- #error "FLT_RADIX != 2. Write your own __printf_fp."
- #endif
-
- #undef alloca /* gmp-impl.h defines it again. */
- #include "gmp.h"
- #include "gmp-impl.h"
- #include "longlong.h"
-
- #ifndef NDEBUG
- static void mpn_dump (const char *str, mp_limb *p, mp_size_t size);
- #define MPN_DUMP(x,y,z) mpn_dump(x,y,z)
- #else
- #define MPN_DUMP(x,y,z)
- #endif
-
- extern mp_size_t __mpn_extract_double (mp_ptr res_ptr, mp_size_t size,
- int *expt, int *is_neg,
- double value);
-
- /* We believe that these variables need as many bits as the largest binary
- exponent of a double. But we are not confident, so we add a few words. */
- #define MPNSIZE ((DBL_MAX_EXP + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) + 3
-
- #define MPN_VAR(name) mp_limb name[MPNSIZE]; mp_size_t name##size
- #define MPN_ASSIGN(dst,src) \
- memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb))
- #define MPN_POW2(dst, power) \
- do { \
- MPN_ZERO (dst, (power) / BITS_PER_MP_LIMB); \
- dst[(power) / BITS_PER_MP_LIMB] = \
- (mp_limb) 1 << (power) % BITS_PER_MP_LIMB; \
- dst##size = (power) / BITS_PER_MP_LIMB + 1; \
- } while (0)
-
- /* Compare *normalized* mpn vars. */
- #define MPN_GT(u,v) \
- (u##size > v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) > 0))
- #define MPN_LT(u,v) \
- (u##size < v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) < 0))
- #define MPN_GE(u,v) \
- (u##size > v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) >= 0))
- #define MPN_LE(u,v) \
- (u##size < v##size || (u##size == v##size && __mpn_cmp (u, v, u##size) <= 0))
- #define MPN_EQ(u,v) \
- (u##size == v##size && __mpn_cmp (u, v, u##size) == 0)
- #define MPN_NE(u,v) \
- (!MPN_EQ(u,v))
-
- int
- DEFUN(__printf_fp, (s, info, args),
- FILE *s AND CONST struct printf_info *info AND va_list *args)
- {
- mp_limb cy;
-
- int done = 0;
-
- /* Decimal point character. */
- CONST char *CONST decimal = _numeric_info->decimal_point;
-
- LONG_DOUBLE fpnum; /* Input. */
- int is_neg;
-
- MPN_VAR (f); /* Fraction. */
-
- int e; /* Base-2 exponent of the input. */
- CONST int p = DBL_MANT_DIG; /* Internal precision. */
- MPN_VAR (scale); MPN_VAR (scale2); MPN_VAR (scale10); /* Scale factor. */
- MPN_VAR (loerr); MPN_VAR (hierr); /* Potential error in the fraction. */
- int k; /* Digits to the left of the decimal point. */
- int cutoff; /* Where to stop generating digits. */
- MPN_VAR (r); MPN_VAR (r2); MPN_VAR (r10); /* Remainder. */
- int roundup;
- int low, high;
- char digit;
-
- MPN_VAR (tmp); /* Scratch space. */
-
- int j;
-
- char type = tolower (info->spec);
- int prec = info->prec;
- int width = info->width;
-
- /* This algorithm has the nice property of not needing a buffer.
- However, to get the padding right for %g format, we need to know
- the length of the number before printing it. */
-
- #ifndef LDBL_DIG
- #define LDBL_DIG DBL_DIG
- #endif
- #ifndef LDBL_MAX_10_EXP
- #define LDBL_MAX_10_EXP DBL_MAX_10_EXP
- #endif
-
- char *buf = __alloca ((prec > LDBL_DIG ? prec : LDBL_DIG) +
- LDBL_MAX_10_EXP + 3); /* Dot, e, exp. sign. */
- register char *bp = buf;
- #define put(c) *bp++ = (c)
-
-
- /* Produce the next digit in DIGIT.
- Return nonzero if it is the last. */
- inline int hack_digit (void)
- {
- int cnt;
- mp_limb high_qlimb;
-
- --k;
- cy = __mpn_mul_1 (r10, r, rsize, 10);
- r10size = rsize;
- if (cy != 0)
- r10[r10size++] = cy;
-
- MPN_DUMP ("r", r, rsize);
- MPN_DUMP ("r10", r10, r10size);
- MPN_DUMP ("scale", scale, scalesize);
-
- /* Compute tmp = R10 / scale and R10 = R10 % scale. */
- count_leading_zeros (cnt, scale[scalesize - 1]);
- if (cnt != 0)
- {
- mp_limb norm_scale[scalesize];
- mp_limb cy;
- assert (scalesize != 0);
- __mpn_lshift (norm_scale, scale, scalesize, cnt);
- assert (r10size != 0);
- cy = __mpn_lshift (r10, r10, r10size, cnt);
- if (cy != 0)
- r10[r10size++] = cy;
- high_qlimb = __mpn_divmod (tmp, r10, r10size, norm_scale, scalesize);
- tmp[r10size - scalesize] = high_qlimb;
- r10size = scalesize;
- __mpn_rshift (r10, r10, r10size, cnt);
- }
- else
- {
- high_qlimb = __mpn_divmod (tmp, r10, r10size, scale, scalesize);
- tmp[r10size - scalesize] = high_qlimb;
- r10size = scalesize;
- }
-
- MPN_DUMP ("high_qlimb", &high_qlimb, 1);
- MPN_DUMP ("r10", r10, r10size);
-
- /* We should have a quotient < 10. It might be stored */
- high_qlimb = tmp[0];
- digit = '0' + high_qlimb;
-
- r10size = __mpn_normal_size (r10, r10size);
- if (r10size == 0)
- /* We are not prepared for an mpn variable with zero limbs. */
- r10size = 1;
-
- MPN_ASSIGN (r, r10);
- assert (rsize != 0);
- cy = __mpn_lshift (r2, r, rsize, 1);
- r2size = rsize;
- if (cy != 0)
- r2[r2size++] = cy;
-
- cy = __mpn_mul_1 (loerr, loerr, loerrsize, 10);
- if (cy)
- loerr[loerrsize++] = cy;
- cy = __mpn_mul_1 (hierr, hierr, hierrsize, 10);
- if (cy)
- hierr[hierrsize++] = cy;
-
- low = MPN_LT (r2, loerr);
-
- /* tmp = scale2 - hierr; */
- if (scale2size < hierrsize)
- high = 1;
- else
- {
- cy = __mpn_sub (tmp, scale2, scale2size, hierr, hierrsize);
- tmpsize = scale2size;
- high = cy || (roundup ? MPN_GE (r2, tmp) : MPN_GT (r2, tmp));
- }
-
- if (low || high || k == cutoff)
- {
- /* This is confusing, since the text and the code in Steele's and
- White's paper are contradictory. Problem numbers:
- printf("%20.15e\n", <1/2^106>) is printed as
- 1.232595164407830e-32 (instead of 1.232595164407831e-32)
- if we obey the description in the text;
- 1/2^330 is badly misprinted if we obey the code. */
- if (high && !low)
- ++digit;
- #define OBEY_TEXT 1
- #if OBEY_TEXT
- else if (high && low && MPN_GT (r2, scale))
- #else
- else if (high == low && MPN_GT (r2, scale))
- #endif
- ++digit;
- return 1;
- }
-
- return 0;
- }
-
- const char *special = NULL; /* "NaN" or "Inf" for the special cases. */
-
- /* Fetch the argument value. */
- if (info->is_long_double)
- fpnum = va_arg (*args, LONG_DOUBLE);
- else
- fpnum = (LONG_DOUBLE) va_arg (*args, double);
-
- /* Check for special values: not a number or infinity. */
-
- if (__isnan ((double) fpnum))
- {
- special = "NaN";
- is_neg = 0;
- }
- else if (__isinf ((double) fpnum))
- {
- special = "Inf";
- is_neg = fpnum < 0;
- }
-
- if (special)
- {
- int width = info->prec > info->width ? info->prec : info->width;
-
- if (is_neg || info->showsign || info->space)
- --width;
- width -= 3;
-
- if (!info->left)
- while (width-- > 0)
- outchar (' ');
-
- if (is_neg)
- outchar ('-');
- else if (info->showsign)
- outchar ('+');
- else if (info->space)
- outchar (' ');
-
- {
- register size_t len = 3;
- while (len-- > 0)
- outchar (*special++);
- }
-
- if (info->left)
- while (width-- > 0)
- outchar (' ');
-
- return done;
- }
-
- /* Split the number into a fraction and base-2 exponent. The fractional
- part is scaled by the highest possible number of significant bits of
- fraction. We represent the fractional part as a (very) large integer. */
-
- fsize = __mpn_extract_double (f, sizeof (f) / sizeof (f[0]),
- &e, &is_neg, fpnum);
-
- if (prec == -1)
- prec = 6;
- else if (prec == 0 && type == 'g')
- prec = 1;
-
- if (type == 'g')
- {
- if (fpnum != 0)
- {
- if (is_neg)
- fpnum = - fpnum;
-
- if (fpnum < 1e-4)
- type = 'e';
- else
- { /* XXX do this more efficiently */
- /* Is (int) floor (log10 (FPNUM)) >= PREC? */
- LONG_DOUBLE power = 10;
- j = prec;
- if (j > p)
- j = p;
- while (--j > 0)
- {
- power *= 10;
- if (fpnum < power)
- /* log10 (POWER) == floor (log10 (FPNUM)) + 1
- log10 (FPNUM) == Number of iterations minus one. */
- break;
- }
- if (j <= 0)
- /* We got all the way through the loop and F (i.e., 10**J)
- never reached FPNUM, so we want to use %e format. */
- type = 'e';
- }
- }
-
- /* For 'g'/'G' format, the precision specifies "significant digits",
- not digits to come after the decimal point. */
- --prec;
- }
-
- if (fsize == 1 && f[0] == 0)
- /* Special case for zero.
- The general algorithm does not work for zero. */
- {
- put ('0');
- if (tolower (info->spec) != 'g' || info->alt)
- {
- if (prec > 0 || info->alt)
- put (*decimal);
- while (prec-- > 0)
- put ('0');
- }
- if (type == 'e')
- {
- put (info->spec);
- put ('+');
- put ('0');
- put ('0');
- }
- }
- else
- {
- cutoff = -prec;
-
- roundup = 0;
-
- if (e > p)
- {
- /* The exponent is bigger than the number of fractional digits. */
- MPN_ZERO (r, (e - p) / BITS_PER_MP_LIMB);
- if ((e - p) % BITS_PER_MP_LIMB == 0)
- {
- MPN_COPY (r + (e - p) / BITS_PER_MP_LIMB, f, fsize);
- rsize = fsize + (e - p) / BITS_PER_MP_LIMB;
- assert (rsize != 0);
- }
- else
- {
- assert (fsize != 0);
- cy = __mpn_lshift (r + (e - p) / BITS_PER_MP_LIMB, f, fsize,
- (e - p) % BITS_PER_MP_LIMB);
- rsize = fsize + (e - p) / BITS_PER_MP_LIMB;
- if (cy)
- r[rsize++] = cy;
- }
-
- MPN_POW2 (scale, 0);
- assert (scalesize != 0);
-
- /* The number is (E - P) factors of two larger than
- the fraction can represent; this is the potential error. */
- MPN_POW2 (loerr, e - p);
- assert (loerrsize != 0);
- }
- else
- {
- /* The number of fractional digits is greater than the exponent.
- Scale by the difference factors of two. */
- MPN_ASSIGN (r, f);
- MPN_POW2 (scale, p - e);
- MPN_POW2 (loerr, 0);
- }
- MPN_ASSIGN (hierr, loerr);
-
- /* Fixup. */
-
- MPN_POW2 (tmp, p - 1);
-
- if (MPN_EQ (f, tmp))
- {
- /* Account for unequal gaps. */
- assert (hierrsize != 0);
- cy = __mpn_lshift (hierr, hierr, hierrsize, 1);
- if (cy)
- hierr[hierrsize++] = cy;
-
- assert (rsize != 0);
- cy = __mpn_lshift (r, r, rsize, 1);
- if (cy)
- r[rsize++] = cy;
-
- assert (scalesize != 0);
- cy = __mpn_lshift (scale, scale, scalesize, 1);
- if (cy)
- scale[scalesize++] = cy;
- }
-
- /* scale10 = ceil (scale / 10.0). */
- if (__mpn_divmod_1 (scale10, scale, scalesize, 10) != 0)
- {
- /* We got a remainder. __mpn_divmod_1 has floor'ed the quotient
- but we want it to be ceil'ed. Adjust. */
- cy = __mpn_add_1 (scale10, scale10, scalesize, 1);
- if (cy)
- abort ();
- }
- scale10size = scalesize;
- scale10size -= scale10[scale10size - 1] == 0;
-
- k = 0;
- while (MPN_LT (r, scale10))
- {
- mp_limb cy;
-
- --k;
-
- cy = __mpn_mul_1 (r, r, rsize, 10);
- if (cy != 0)
- r[rsize++] = cy;
-
- cy = __mpn_mul_1 (loerr, loerr, loerrsize, 10);
- if (cy != 0)
- loerr[loerrsize++] = cy;
-
- cy = __mpn_mul_1 (hierr, hierr, hierrsize, 10);
- if (cy != 0)
- hierr[hierrsize++] = cy;
- }
-
- do
- {
- mp_limb cy;
- assert (rsize != 0);
- cy = __mpn_lshift (r2, r, rsize, 1);
- r2size = rsize;
- if (cy != 0)
- r2[r2size++] = cy;
-
- /* tmp = r2 + hierr; */
- if (r2size > hierrsize)
- {
- cy = __mpn_add (tmp, r2, r2size, hierr, hierrsize);
- tmpsize = r2size;
- }
- else
- {
- cy = __mpn_add (tmp, hierr, hierrsize, r2, r2size);
- tmpsize = hierrsize;
- }
- if (cy != 0)
- tmp[tmpsize++] = cy;
-
- /* while (r2 + hierr >= 2 * scale) */
- assert (scalesize != 0);
- cy = __mpn_lshift (scale2, scale, scalesize, 1);
- scale2size = scalesize;
- if (cy)
- scale2[scale2size++] = cy;
- while (MPN_GE (tmp, scale2))
- {
- cy = __mpn_mul_1 (scale, scale, scalesize, 10);
- if (cy)
- scale[scalesize++] = cy;
- ++k;
- assert (scalesize != 0);
- cy = __mpn_lshift (scale2, scale, scalesize, 1);
- scale2size = scalesize;
- if (cy)
- scale2[scale2size++] = cy;
- }
-
- /* Perform any necessary adjustment of loerr and hierr to
- take into account the formatting requirements. */
-
- if (type == 'e')
- cutoff += k - 1; /* CutOffMode == "relative". */
- /* Otherwise CutOffMode == "absolute". */
-
- { /* CutOffAdjust. */
- int a = cutoff - k;
- MPN_VAR (y);
- MPN_ASSIGN (y, scale);
-
- /* There is probably a better way to do this. */
-
- while (a > 0)
- {
- cy = __mpn_mul_1 (y, y, ysize, 10);
- if (cy)
- y[ysize++] = cy;
- --a;
- }
- while (a < 0)
- {
- if (__mpn_divmod_1 (y, y, ysize, 10) != 0)
- {
- /* We got a remainder. __mpn_divmod_1 has floor'ed the
- quotient but we want it to be ceil'ed. Adjust. */
- cy = __mpn_add_1 (y, y, ysize, 1);
- if (cy)
- abort ();
- }
- ysize -= y[ysize - 1] == 0;
- ++a;
- }
-
- if (MPN_GT (y, loerr))
- MPN_ASSIGN (loerr, y);
- if (MPN_GE (y, hierr))
- {
- MPN_ASSIGN (hierr, y);
- roundup = 1;
- /* Recalculate: tmp = r2 + hierr */
- if (r2size > hierrsize)
- {
- cy = __mpn_add (tmp, r2, r2size, hierr, hierrsize);
- tmpsize = r2size;
- }
- else
- {
- cy = __mpn_add (tmp, hierr, hierrsize, r2, r2size);
- tmpsize = hierrsize;
- }
- if (cy != 0)
- tmp[tmpsize++] = cy;
- }
- } /* End CutOffAdjust. */
-
- } while (MPN_GE (tmp, scale2));
-
- /* End Fixup. */
-
- /* First digit. */
-
- hack_digit ();
-
- if (type == 'e')
- {
- /* Exponential notation. */
-
- int expt = k; /* Base-10 exponent. */
- int expt_neg;
-
- expt_neg = k < 0;
- if (expt_neg)
- expt = - expt;
-
- /* Find the magnitude of the exponent. */
- j = 10;
- while (j <= expt)
- j *= 10;
-
- /* Write the first digit. */
- put (digit);
-
- if (low || high || k == cutoff)
- {
- if ((tolower (info->spec) != 'g' && prec > 0) || info->alt)
- put (*decimal);
- }
- else
- {
- int stop;
-
- put (*decimal);
-
- /* Remaining digits. */
- do
- {
- stop = hack_digit ();
- put (digit);
- } while (! stop);
- }
-
- if (tolower (info->spec) != 'g' || info->alt)
- /* Pad with zeros. */
- while (--k >= cutoff)
- put ('0');
-
- /* Write the exponent. */
- put (isupper (info->spec) ? 'E' : 'e');
- put (expt_neg ? '-' : '+');
- if (expt < 10)
- /* Exponent always has at least two digits. */
- put ('0');
- do
- {
- j /= 10;
- put ('0' + (expt / j));
- expt %= j;
- }
- while (j > 1);
- }
- else
- {
- /* Decimal fraction notation. */
-
- if (k < 0)
- {
- put ('0');
- if (prec > 0 || info->alt)
- put (*decimal);
-
- /* Write leading fractional zeros. */
- j = 0;
- while (--j > k)
- put ('0');
- }
-
- put (digit);
- if (!low && !high && k != cutoff)
- {
- int stop;
- do
- {
- stop = hack_digit ();
- if (k == -1)
- put (*decimal);
- put (digit);
- } while (! stop);
- }
-
- while (k > 0)
- {
- put ('0');
- --k;
- }
- if ((type != 'g' && prec > 0) || info->alt)
- {
- if (k == 0)
- put (*decimal);
- while (k-- > -prec)
- put ('0');
- }
- }
- }
-
- #undef put
-
- /* The number is all converted in BUF.
- Now write it with sign and appropriate padding. */
-
- if (is_neg || info->showsign || info->space)
- --width;
-
- width -= bp - buf;
-
- if (!info->left && info->pad == ' ')
- /* Pad with spaces on the left. */
- while (width-- > 0)
- outchar (' ');
-
- /* Write the sign. */
- if (is_neg)
- outchar ('-');
- else if (info->showsign)
- outchar ('+');
- else if (info->space)
- outchar (' ');
-
- if (!info->left && info->pad == '0')
- /* Pad with zeros on the left. */
- while (width-- > 0)
- outchar ('0');
-
- if (fwrite (buf, bp - buf, 1, s) != 1)
- return -1;
- done += bp - buf;
-
- if (info->left)
- /* Pad with spaces on the right. */
- while (width-- > 0)
- outchar (' ');
-
- return done;
- }
-
- #ifndef NDEBUG
- static void
- mpn_dump (str, p, size)
- const char *str;
- mp_limb *p;
- mp_size_t size;
- {
- fprintf (stderr, "%s = ", str);
- while (size != 0)
- {
- size--;
- fprintf (stderr, "%08lX", p[size]);
- }
- fprintf (stderr, "\n");
- }
- #endif
-